Constraint Molecular Dynamics approach to Fermionic systems. 
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We propose a Constraint Molecular Dynamics model for Fermionic system. In this approach 
the equations of motion of wave packets for the nuclear many-body problem are solved by imposing 
that the one-body occupation probability f(r,p,t) can assume only values less or equal to 1. This 
condition reflects the Fermionic nature of the studied systems and it is implemented with a fast 
algorithm which allows also the study of the heaviest colliding system. The parameters of the model 
have been chosen to reproduce the average binding energy and radii of nuclei in the mass region 
A = 30 ~ 208. Some comparison to data is given. 
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D . I. INTRODUCTION 

Q 

. Heavy ion collisions in the medium energy region have been described in a large variety of semiclassical approaches 
' to the many-body problem. As it is well known the one-body semiclassical transport models like the Boltzman- 
Northeim-Vlasov (BNV) @ and Vlasov-Uehling-Uhlenbeck (VUU) [§ are not suited to describe processes in which 
t-H ■ a large number of final fragments arc produced. This is due to the fact that the correlations treated in the one-body 
^ | approach are not able to describe the large fluctuations which develops in a multifragmentation process. 

This difficulty can be solved by adopting more suitable treatments of the A-body problem like molecular dynamics. 
Several molecular dynamics models have been developed up to now jlj. In the quantum molecular dynamics (QMD) jij] 
the A-body wave function is expressed through a direct product of wave packets each of which satisfies the minimum 
1 uncertainty relation cr r a p = h/2, where o r and a p represent the dispersion of the corresponding Wigner transform in 
f^i , configuration and momentum space respectively. 

The Fermionic nature of the nuclear many-body problem has been considered in the Fermionic molecular dynamics 
H and more recently in the antisymmetrized molecular dynamics (AMD) || . In these models, the wave function of the 
• system is expressed as a single Slater determinant of A wave packets. In this way the Fermionic nature of the system is 
preserved. In particular in the AMD approach two-body collisions are introduced in the "physical coordinates" which 
are obtained by a canonical transformation from the parameter coordinates of wave-packets Q . The nucleon-nucleon 
collision is only allowed when an inverse transformation from a newly chosen "physical coordinate" (the final state of 
the collision) into the corresponding parameter coordinate exists. This condition is understood as a stochastic change 
of the state from a single Slater determinant into another, both of which have occupation probability 1 or 0. Due to 
the four dimensional matrix element of two-body interaction, the CPU time necessary to work out calculations for 
5J! ■ systems with total mass larger than 200 is very large for practical studies. 

On the other hand the QMD calculations are much easier to carry out because they need in general only double-fold 
loops to calculate two-body interactions. But it is obvious that the Fermionic feature is lacking in the QMD model. 
Two-body collisions are also introduced phenomenologically. The Pauli blocking of the final state is usually checked 
in a similar manner to the collision term in the BNV model. The two-body collisions with Pauli blocking have some 
effects to maintain the Fermionic feature of the system. However, in the ground states or in low energy reaction 
phenomena, two-body collisions are absent or very rare. Even if the initial state is in good agreement with the phase 
space distribution of a Fermionic system, the time evolution by the classical equations of motion surely breaks the 
initial distribution which evolves to a classical Boltzmann one. 

To compensate this shortcoming, the Pauli potential is introduced by several authors to mimic the Fermionic 
features [^|-^| • This phenomenological potential forbids nucleons of the same spin and isospin from coming close to 
each other in the phase space. Although the Pauli potential gives us some good results such as stable ground states, 
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with energies in good agreement with experiments and saturation properties of nuclear matter, it also gave some 
undesirable by-products, for instance, spurious repulsion in the collision problem. 

In the present work we propose a new molecular dynamics model: the Constraint Molecular Dynamics (CoMD), 
which aims to overcome the above mentioned limitations. In particular we want: 

i) to describe the Fermionic nature of the TV-body system with the more general condition that the occupation 
probability / < 1; 

ii) to realize a model for which the computational time is short enough to allow the study of the heaviest systems. 



II. THE MODEL 



A. Theoretical framework 



In QMD, each nucleon state is represented by a Gaussian wave- function of width ay, 
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where (r^) and (pi) are the centers of position and momentum of i-th nucleon, respectively. The total wave-function 
is assumed to be a direct product of these wave- functions. Therefore the N-body distribution function is the direct 
product of the single particle distribution functions fi. The fi are obtained by the Wigner transform of the wave- 
functions 4>i and the one-body distribution function is given by: 
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We note that since ay is a real number in the QMD approach the distribution function fi produce the minimum 
uncertainty relation a r a. p = h/2 in one-body phase space. 

In this paper we extend this kind of representation. We in fact take the dispersion in momentum a p as a parameter 
as well as that in coordinate space allowing for the more general uncertainty condition a p a r > h/2. Therefore we 
write the one-body distribution function as: 
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This extension implies that each single particle wave function fa can be now represented by an incoherent super- 
imposition of a large number of minimum uncertainty wave-packets. 

The equations of motion of (rj) and (pj) arc derived using the time-dependent variational principle which gives: 
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The Hamiltonian H consists of the kinetic energy and the two-body effective interactions: 

(P,) 2 



Va = 



2m 2 

;3 



= J d 3 r t tfrj PiiTijpjiTjJVij, 



Pr 



d 3 p /i(r,p), 



(5) 

(6) 

(7) 
(8) 



where Vij represents two body interaction between particles i and j and pi represents the density distribution of 
nucleon i. One should note that the kinetic energy coming from the momentum dispersion 3a p 2 /2m is omitted as a 
spurious constant term. 
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B. Effective interaction 
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The nucleon-nucleon interaction has been described through the sum of several terms: 

Vij = V ml + y (3) + V asy + V sup + V Coul (9) 

(10) 

(11) 

L] (12) 

(13) 
(14) 

In the above relations the power operator [..] 7 ^ 6 and V 2 act after the folding operations indicated in Eq. (Q). The 
Tj indicates the isospin degree of freedom. The V vo1 and the terms represent the two-body and the so called 
three-body term of the well known Skyrme interaction. The value of to and £3 have been fixed to —356 MeV and 
303 MeV. These values reproduce the saturation density po and binding energy for symmetric nuclear matter with a 
compressibility of 200 MeV. The third term represent the asymmetry term with a sym = 32 MeV. The parameter C s 
of the "surface" term V surf is chosen through a global fit procedure to reproduce, on the average, the radii and the 
binding energies for nuclei with mass-ranging from 30 to 208. 
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C. Numerical methods and the Constraint 

The set of equations expressed in Eq. (^|) have been solved using a fourth-order Runge-Kutta method coupled with 
two numerical algorithms which aim to take into account the effects of the residual interaction and the Fermionic 
nature of the many-body problem we are studying. 

One consists of the usual two-body elastic collisions which mimics the effect of the short range repulsive residual 
interaction together with the stochastic change of the phase-space distribution with Pauli blocking in the final states. 
The isospin dependent parameterization of the nucleon-nucleon elastic angular distribution together with the concept 
of mean-free path have been used to compute the collision probability per unit of time Q. The Pauli blocking factor, 
which is related to the constraint, is discussed later. 

The other algorithm constraints at each time step the following quantities: 

7i < 1 (for all i) (15) 

7i = E<W^ / /;( r -p; < r A (p.)) d * r d3 p (i 6 ) 

The coordinate Si represent the nucleon spin projection quantum number. The integral is performed on an hypercube 
of volume h 3 in the phase-space centered around the point ((i"0, (p;)) with size J cr r and y^^-Op in the r and 

p spaces respectively. The quantities / i can therefore be interpreted like an occupation probability of the one-body 
phase-space around the point of coordinate ((r^), (p^)). In general more than 60% of the occupation probability f i 
arises from the contribution fi of the particle i itself. To realize the constraint expressed trough the relation (|l^) the 
following procedure has been used: 

At each time step and for each particle i an ensemble Ki of nearest identical particles (including the particle i) is 
determined within the distances 3a r and 3a p in the phase space. If the phase space occupation f i has a value greater 
than 1 we perform many body elastic scattering among the ensemble Ki in order to reduce the phase space occupation. 
To reduce the CPU time the many body scattering process is practically done by a series of two-body scatterings. 
The Pauli blocking probability for the usual two-body collision process is given as: 

i\,iock = l- (1-50(1-5,-), (17) 
S ; e min I 1, 1 . ( LS) 
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fi = / .A:(r,p;(r 4 ),( Pl }) d 3 r d 3 p. 
Jh 3 



(19) 



The quantity Si in Eq. ([18|) is the "effective occupation" of phase space at the position of particle i. In fact in the 
definition of Si we are subtracting the contribution of the particle i itself. 

We stress that the constraint acts in a way complementary to the collision term. In fact particles of low momenta 
are strongly effected by the constraint in such a way to avoid that the distribution becomes a classical one. On the 
contrary the collision term is especially important for particles located at high relative momenta. 



III. CALCULATIONS AND COMPARISON WITH OTHER MODELS 



As an example we compare our results on the isotope distribution in 40 Ca + 40 Ca at 35 MeV/nucleon with the 
experime nt Eofl and to the result of QMD. Finally we compare CoMD results to experimental data on central Au+Au 
collisions [|11||. 



A. Initialization 



The ground state configuration of the nuclei is obtained using a modified cooling procedure. The nucleon are first 
distributed in a sphere of radius 1.4 x A 1 / 3 fm in configuration space and in a sphere of radius Pp m (Fermi momentum 
for infinite nuclear matter) in momentum space. The equations of motion with friction terms are solved coupled with 
the constraint. At each time step if the value of f i is greater than 1 the momenta of the particles belonging to the 
ensemble Ki are scaled by a factor 1.02. If f i is less than 1 the scale factor is set to 0.98. With this procedure the 
total energy and the radius R of the nuclei will reach some stationary values. If the deviations of these values from 
the experimental ones are within 7%, a check is done on the stability. At this stage the friction term for the "cooling" 
is switched off and the time dependence of the nuclear radius R is checked. If R is stable at least for 1000 fm/c 
(inside the confidence level), this initial condition is accepted. In Fig. 1 we show R as function of time for two typical 
"good events" representing the 40 Ca and 208 Pb nuclei. The binding energies are —8.2 and —8.4 MeV respectively. 
The reason of the stability in the CoMD case is the constraint. 
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FIG. 1. Nuclear radii R as function of time for typical ground state configurations of 40 Ca and 208 Pb nuclei. 

In fact without constraint, the initial Fcrmionic distribution will evolve to a classical one. This means that during 
the interaction several particles will tend to have low relative momenta trying to reach a stable classical configuration 
(a classical ground state would be a solid) while some other particles (as due to the energy conservation) will acquire 
momenta large enough to overcome the barrier and leave the system. This does not happen in CoMD since the 
constraint maintains the Fermionic nature of the system. 

The results shown in Fig. 1 are obtained by setting ay — 1.3 fm, a p /h — 0.47 fm _1 and C s = —2 MeV-fm 2 . 
The value of the surface term might be surprising at first, but we must notice that part of the surface term comes 
automatically when using Gaussian. We observe also that using this procedure with the above set of parameters the 
effective average kinetic energy per nucleon K (the kinetic energy regarding the centroid of the Gaussian wave packet) 
is very close to the value obtained in a Thomas Fermi model (about 20 MeV) This feature is instead lost in most of 
the molecular dynamics models while it is a peculiarity of the proposed one. 

In fact in the present approach for a fixed average density and effective interaction, the ratio between kinetic energy 
and the potential one crucially depends on a p which can be varied to values greater than h/2a r (see Sec. II). In 
particular a v determines the filling of the phase-space , and therefore, the constrained value of / will be automatically 
reached with a value of K strongly depending on a p . 

This feature will allow us, in the global fit procedure, to find a set of parameters which reproduce the nuclear 
binding energies, the radii and the ratio between K and the potential energy close to the value of the Thomas Fermi 
model. 
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B. The Nucleus-Nucleus Collision 



In this section we show some results concerning the 40 Ca + 40 Ca collision at £i a b = 35 MeV/nucleon. In Fig. 2 
we display the time dependence of f i for all particle index i in a typical event with zero impact parameter. The 
upper and the bottom parts show the results of QMD and CoMD models, respectively, with exactly the same initial 
condition. The difference between the two cases is quite clear. In the QMD case after a short time (about 50 fm/c) 
the values of the f i for most of the particles start to be greater than 1 and are larger than 2 after some hundreds of 
fm/c. This result is easily understood. In the QMD case, which has no constraint regarding the Fcrmionic nature 
except for the two-nucleon collision process, the momentum distribution tends to the classical limit. Obviously in the 
CoMD case, due to the constraint, the phase space occupation f i is on the average less than 1 at each time step. We 
note a good uniformity of j i as function of the particle index i and time t. These results obviously affect the collision 
rate r c . 



QMD 




FIG. 2. Two-dimensional histograms showing the occupation probability fi as function of the particle index i and time for 
the system under study. The upper panel refers to the case QMD case while the lower one is relative to the CoMD model. 

In Fig. 3 we show r c as function of t for central collision for both the QMD and the CoMD. In the QMD case the 
collision rate steeply increases and decreases around 40 fm/c. This behavior can be explained as follows: although we 
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use the same initial condition as CoMD, there appear some fluctuation of occupation number f i due to the classical 
nature of QMD. Therefore two-nucleon collisions are easily to occur at the beginning compared to CoMD. After 
about 30 fm/c, however, when the two nuclei start to overlap the j i increases spuriously above 1. This causes the 
rapid decrease of the collision rate in QMD calculations because of the Pauli blocking. The behavior is completely 
different in the CoMD case. The collision rate develops in about 200 fm/c and reaches a value 3 times greater than 
the maximum value relative to the QMD case, which in turn gives more stopping. After this time the disassembled 
system, i.e. the fragments, gradually approach their ground state, which makes the decrease of the collision rate. 
One may notice, however, the unusually large collision rate even after several hundreds of fm/c. The reason of this 
large number is due to collisions with very low relative momentum which we do not exclude in the calculations. If 
some low-momentum cut off of two-nucleon collision is put, we could avoid those high collision rate at later times. 
Nevertheless those collisions with low relative momentum have no effect on the dynamics. 
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FIG. 3. Collision rate r c as function of time for the central collision 40 Ca + 40 Ca at 35 MeV/nucleon in the QMD and 
CoMD cases. 

In Figs. 4 and 5 we compare the isotope distribution given by QMD and CoMD approaches with the experimental 
data on 40 Ca + 40 Ca |l(]]. About 2000 events have been generated for an impact parameter range 6:0^8 fm. 
The minimum spanning tree method jj| in the configuration space has been applied to determinate the fragments 
at t = 300 fm/c and at t = 3000 fm/c. We have verified that within the time 3000 fm/c all the fragments are 
stable. What we can also see from the figure is that the main features of the fragment mass distribution are already 
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determined at 300 fm/c. After that the process is dominated by the emission of particles from the heavier fragments. 
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FIG. 4. Comparison between the experimental isotope distribution measured for the 40 Ca + 40 Ca at 35 MeV/nucleon [10] 
system and the theoretical prediction performed according to the QMD approach. The calculations are shown at two different 
time intervals as indicated in the figure. 
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FIG. 5. Same as Fig.4 but for CoMD model. 

In the upper panel of Fig. 4 we compare the experimental isotope distribution to QMD at 3000 fm/c. We observe 
a clear disagreement between the QMD calculation and the experiment. The calculation shows essentially binary 
character of the reaction while the experimental data shows a significant production of the intermediate mass fragment 
(IMF: Z > 3). The binary behavior is due to some kind of transparency obtained in the calculations. This in turn is 
due to lack of two-body collisions because of the over crowding of the phase space as discussed above. At the bottom 
of Fig. 4 the result at 300 fm/c is shown. In this case we observe a similar behavior as above, the only difference being 
a shift of the main bump towards higher value of Z and a more pronounced "U" shape. In Fig. 5 the same comparison 
is shown for the CoMD calculation. The calculations are now in a satisfactory agreement with the experimental data. 
In particular we note that the theoretical prediction of the Z = 2 yield is about a factor ten higher than the QMD 
case. This is clearly an effect brought by the constraint (|l5|) which favors a particle states. However the predicted 
yields of Z — 1 and Z = 2 isotope still show a marked difference with the experimental one. 

We have also performed some calculations for the Au+Au system at the same beam energy. Here we are especially 
interested in central collision. In fact it is quite puzzling the behavior of such a small system but so heavily charged. 
In Fig. 6 the charge distribution is given for impact parameters up to 3.5 fm and compared to data JlT[ . A few features 
are worth noticing. 

1) CoMD gives too many protons (not displayed in the figure-out of scale). 

2) The theoretical distribution is slightly shallower than data. Notice however that the data shows a jump for 
Z = 20 due to the detectors used. If we slightly shift the data yield up for Z > 20 we get a better agreement to the 
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CoMD results. 
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FIG. 6. Comparison between the experimental isotope distribution (full circles) [11] and the calculated one according to the 
CoMD model (open circles) for central collisions (b ~ ~ 3.5/m) 197 Au+ 197 Au at 35 MeV/nucleon. The vertical bars indicate 
the errors as due to the statistics counting. 

The system was followed up to 1500 fm/c which should be long enough to get cold fragments. Nevertheless we can 
see from the comparison that the model is working not too bad especially if one recalls that other dynamical model 
calculations give very steep distributions. 

IV. CONCLUSION 

The CoMD model proposed in the present work is able to reproduce with the same set of parameters both the 
main characteristic of stable nuclei in a wide region of mass (A — 30 ~ 208) and the experimental isotope distribution 
produced in the collision 40 Ca + 40 Ca and Au+Au at £iab = 35 MeV/nucleon. This possibility has a considerable 
relevance because the above mentioned inclusive information can not be reproduced by the QMD. To this respect the 
success of the CoMD is due to the constraint represented by the relation (|l5|). This constraint, introduced to describe 
the Fermionic nature of the nuclear many-body problem, affects the dynamics of the nucleus-nucleus collision for two 
main reasons: 

a) the nucleon-nucleon collision rate is higher with respect to the QMD case; 

b) the constraint for low momentum particles produces obviously on the average a non local repulsion effect. 
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Both these effects play a determinant role for the disassembly of the highly excited intermediate systems formed at 
the beginning of a nuclear collision like that investigated in this work. Finally we stress that for the model proposed 
the typical CPU time needed to follow the time evolution of systems of mass number around 80 for 300 fm/c is quite 
short: about 10 sec in a 600 MHz Unix machine. 
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